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We present a formulation of non-Markovian quantum trajectories for open systems from a mea- 
surement theory perspective. In our treatment there are three distinct ways in which non-Markovian 
behavior can arise; a mode dependent coupling between bath (reservoir) and system, a dispersive 
bath, and by spectral detection of the output into the bath. In the first two cases the non-Markovian 
behavior is intrinsic to the interaction, in the third case the non-Markovian behavior arises from 
the method of detection. We focus in detail on the trajectories which simulate real-time spectral 
detection of the light emitted from a localized system. In this case, the non-Markovian behavior 
arises from the uncertainty in the time of emission of particles that are later detected. The results of 
computer simulations of the spectral detection of the spontaneous emission from a strongly driven 
two-level atom are presented. 



Quantum trajectories or stochastic Schrodinger equations |l|-[3| have been used as an effective numerical tool (for 
example |i|-|7j]) to solve the time evolution of the reduced density matrix for an open system. In addition, a single 
quantum trajectory simulates the evolution of a system undergoing continuous measurement of its environment 
P,p[-pXj| . This is of fundamental importance in understanding open system behavior as features of a single realization 
of a measurement process can be obscured by density matrix methods which average over the individual realizations 
(see for a dramatic example). 

Traditionally, quantum trajectories have described Markov processes, where the conditioned state of the system at a 
certain time contains all the information needed to calculate the measurement probabilities over the next infinitesimal 
interval. In this case, a master equation can describe the evolution of the reduced density matrix for the open system. 
Recently, however, there has been interest in studying open systems for which the Markov approximation cannot be 
made. For example, in a photonic band-gap material the correlation time of the electromagnetic vacuum fluctuations 
near a band edge is not small on the timescale of the emitting atomic system p2[ , p!3[ . As another example, a Bose- 
Einstein condensate losing atoms due to a radio-frequency output coupler [ fl4[ experiences Non-Markovian decay for 
certain coupling rates, as the dispersion relation for atoms in free space leads to a non- vanishing correlation time for 
the atoms coupled out of the condensate Jl5| , ^6[ . There are also measurement processes that lead to a non-Markovian 
evolution of the conditioned system p"7| , |18| ]. Real-time spectral measurement is an example Jl9| . In this case, the 
measurement process yields not only temporal information but also frequency information about the system p0| , pl| . 
It is therefore of interest to generalize quantum trajectory techniques to deal with these non-Markovian processes 
p3|-po| . In the non-Markovian case the future measurement probabilities need to be calculated from the evolution of 
the system in the past. 

In this paper we formally derive in Sec.| the evolution equation for a non-Markovian trajectory for arbitrary 
measurement schemes. As an illustrative example we investigate in Sec.|| the trajectories that arise from spectral 
detection of the light from a strongly driven two-level atom. 

A derivation of non-Markovian quantum trajectories from a microscopic model (in the case of heterodyne detection 
p2|) has been made recently by Diosi and Strutz |24|] via path integrals. By rewriting the path integral evolution 
in terms of operators Diosi et al. ]25| ] have used this formulation to simulate trajectories for some simple systems 
with the requirement that the trajectories are consistent with the general formulation. It is not clear (at least to 
the present authors) what physical situations these trajectories correspond to. In contrast, our formulation is from 
a measurement theory perspective. We consider a microscopic model and general measurement schemes that are 
in principle physically realizable. The trajectories are defined from this microscopic model and have an obvious 
interpretation in terms of the measurements. 

A non-Markovian quantum process is one in which the measurement probabilities in the present depend on the 
evolution of the system in the past. This arises naturally when considering the evolution of an open system conditioned 
on measurements made on its emission at a distant detector. In a photon counting experiment, for example, the 
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probability of detecting a photon at the detector can, in general, be expressed as the probability for the system to 
undergo emissions distributed over times in the past. In the Markov case, the probability of a detection is determined 
by that of a single time of emission. In turn, a measurement result at a particular time does not completely determine 
the state of the system at the same time. In a Markov process our knowledge of the state of the system straight 
after the measurement is complete. In the non-Markov case there is a definite evolution that the system must have 
undergone in the past in order to produce the measurement result, but there is always the possibility that future 
measurement results will require us to revise our knowledge of the evolution of the system over this same time. 
However, a conditioned system state can be defined a finite time in the past if we can assume that any events 
occurring in the system earlier than this time will have a negligible effect on the measurement results in the present. 

It is of interest to be able to simulate measurement processes that yield not only temporal information but also 
frequency information about a system. In other words, real-time spectral measurements, where the emission from 
the system is passed through a 'spectrometer' that splits it into separate spectral components before each component 
is measured. Real-time spectral measurement is relevant to situations where different spectral components of the 
emission from a system are time correlated. A well-known example is the three peaked spectrum of the spontaneous 
emission from a strongly driven two- level ato m |2^| . In this system the two side peaks are correlated via photon 
bunching and are individually anti-bunched [^0 27 1. The temporal features are not accessible by methods that time 
average the detections to determine the spectrum. Here we will treat this problem from a real-time spectral detection 
perspective. 

There are two complementary approaches to modeling spectral detection. The first and most physical approach 
is to model the spectrometer as a separate physical system that is being driven by the output from the system of 
interest. The whole extended system of system plus spectrometer can then be evolved forward together as a cascaded 
system p8Lp9|. Imamoglu has recently developed methods using the extended system approach for the more 
general case of a mode-dependent coupling to the bath. The second approach, which perhaps coincides more closely 
with a signal processing treatment, is to consider the spectrometer as a black box where the output from the box 
is related to the input by a spectral response. Here we will consider the second approach as it is more general and 
is, we believe, closer to the idea behind quantum trajectories; of simulating a continuous measurement process at a 
distant detector by evolving a local system forward in time. However, this approach means that we must consider 
non-Markovian evolution of the system, since emission times become more uncertain as the frequency information 
becomes more precise. 



I. NON-MARKOVIAN QUANTUM TRAJECTORIES 



In this first section we outline the general theory of non-Markovian quantum trajectori es f rom a measurement 
theory perspective. The form of the system- bat h coupling that we consider is given in Sec. IB and we introduce a 
complete set of measurement channels in Sec.IC to model measurements that are sensitive to only a range of modes 
of the bath. In Sec. ID we discuss the finite memory-time assumption for the system. 



A. Measurements on an open system 

Consider the typical open system situation of a small localized system surrounded by an infinite bath. The Hamil- 
tonian of the bath and system can be written as the sum of three parts 

H = H sys + ifbath + -Hint) (1) 

where H sys is the Hamiltonian of the free system, -ffbath the Hamiltonian of the free bath and -ff; n t the interaction 
between the two. 

An open system can be thought of in terms of inputs and outputs : the initial bath propagates in towards the 
system; they interact; and the bath propagates away again. We would like to model the situation where measurements 
are being made continuously on the output bath. The probability of getting a particular string of measurement results 
R(t,to) during the interval [Mo), denoted here by P[R(t 7 t )], is given in terms of the initial state by 

P[R(t,t )} =Tr{\R(t,t ))(R(t,t )\p(t)} , (2) 
= Tr{\R(t,t )}{R(t,to)\U tat (t,t )p{t )Ul t (t,t )^ , (3) 

where /7i n t(i, to) is the unitary evolution operator (in the interaction picture) determined by the total Hamiltonian, 
Eq. (Ill) . Note that we will often use the labels R and / to label both a general element of the set of all possible values 
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and also to represent a particular element of the set. We assume that the system and bath are initially uncorrelated 
and we can write the initial density matrix in the form p(to) = Pbath(io)Psys(io)- The density matrix of the bath can 
be written pbath(^o) = P[I(t,to)]\I(t,to))(I(t,to)\ where I(t,to) labels a particular pure state of the input bath. 
Similarly, the system density matrix can be written, p sys (ta) = J\ PbPk(to)]\'>pk{to)}(' l Pk(to)\. We can now write the 
probability of the result R(t,to) in the form 



P[R(t, to)] = P[R(t> to)I{t, to)fl>k(to)], (4) 

Ik 

= ]T P[Rfa t )\I(t, t )ip k (to)}P[I(t, t )]P[Mto)], (5) 

Ik 

where we have used Bayes theorem. From Eq.(^) the conditional probability can be written 

P[R(t,t )\I(t,t )] = (yj(t )\rt RI (t,t )n RI (t,t )\i;(t )), (6) 

where we have defined the operation ]3(| on the system space corresponding to these measurement results over the 
interval [t, to) as 

n R i(t,t ) = (R(t,t Q )\U int (t,t Q )\I(t,t Q )). (7) 

Note that we have suppressed the explicit reference to the probability distribution over the initial system states for 
ease of notation. This is equivalent to assuming that the system was initially in a pure state. The explicit form of 



the operation Qri for a linear, particle-conserving coupling will be given in Sec. IB. For now, we simply assume that 
an explicit form exists. The idea behind a quantum trajectory simulation is expressed by Eq.(^|). If in a number of 
runs of the simulation we choose the initial states randomly with a frequency corresponding to the initial probability 
distribution P[I(t,to)] then the conditional probabilities P[R(t, to)\I(t,to)] of each run will automatically generate 
values of R(t,to) consistent with P[R(t, to)]. 

In general, the above operation will not satisfy the semi-group property characteristic of a Markov process, i.e, 
fl(t,to) Q,(t, s)Q,(s,to) for all s (we have dropped the subscripts for clarity). However, the operation can always be 
written in the form 

n(t,to) = n(t,»)ofi(8,to), (8) 

where we have defined the o-product by 

n(t, s) o n( a , t ) = n k (t, s)n k ( s , t ). (9) 

This is because one can always insert a sum over a complete set of states, 1 = ^ at any time in the unitary 

evolution. The summation in Eq.(j^) need not be over the complete set of states, but is over a subset of these that are 
consistent with both the input and resultant states and therefore depends on the time s. The inability to factorizc 
this operation is due to the fact that the system is entangled with the bath at time s, so that a measurement of the 
output at this time does not completely determine the state of the system at the same time. 

To understand this we need to consider in more detail the differences between Markovian and non-Markovian 
evolution of open systems. Firstly, it is not possible to write down a master equation for the reduced system density 
matrix p sys (t) = Trb a th{/o(^)}- This is evident even to first order perturbation in the interaction Hamiltonian during 
the time interval T m , 



p{t) - P (t - T m ) = / ds[H int (s),p(t-T m )}. (10) 

where H lnt (s) is in the interaction picture and we have assumed that the contributions from the state at times before 
t — T m are negligible. We now write the density matrix at time t — T m in terms of a factorized part and an unfactorized 
part p(t — T m ) — p® ys (t — T m )/9bath(£o) + x(t ~ T m ), where x(t — T m ) contains the entanglement. Inserting this into 
the above equation and taking the trace over the bath variables we find 

p sys (t) - p sy s(i - Tm) = Tr ba th ds[H int (s), x(t - Tm)} | , (11) 
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where we have assumed that Trb a th{^int(s)pbath(io)} = 0, as the mean bath interaction energy can always be taken 
inside the Hamiltonian for the system. This demonstrates that (even to first order in the perturbation) the change 
in the density matrix of the system at time t is determined by the correlations that existed between the bath and 
the system at t — T m . Note that the master equation for the system for a Markov process is derived by going to 
second-order in the perturbation (correlations can only arise from an initial factorized density matrix by going to at 
least second order in the perturbation) and assuming that the correlations decay instantaneously T m — > on the time 
scale of the inverse coupling strength (see for example |u]]). 

Another consequence of the above mentioned entanglement is that to determine the probabilities for a particular 
measurement result at time t we need to consider the evolution of the system from a time t — T m in the past up until 
time t. In addition, even if the initial state of the system is pure it will not stay pure under evolution conditioned on 
continuous measurement of the output, unlike the Markov case ||||. 

Consider the simplest case of continuously counting the emitted particles (e.g., a photodetector counting photons). 
Suppose that a particle that is detected at time t could only have been emitted during the time interval [t, t — T m ), with 
a specific weighting for each possible emission time. Therefore, to determine the probability of a detection at time t 
we need to propagate the system undergoing emissions over the time interval [t, t — T m ). A diagram of this situation is 
given in Fig.^, showing a possible weighting distribution for each emission time for a particle detected at time t. Now 
assume that we have a measurement record over the interval [t,t — T m ) where a single particle was detected at time 
t > t\ > t — T m , (see Fig{l]). We have the complete measurement record from t to t — T m . Therefore, the conditioned 
state of the emitting system at time t — T m cannot change due to future detections. We are interested, then, in the 
conditioned state at time t — T m . The detection at time t\ corresponds to two possibilities at time t — T m , either the 
particle was emitted after time t — T m or before. In the Hilbert space of the system and bath this corresponds to a 
superposition of the two possibilities: a particle has been as emitted and is present in the bath and a particle has not 
been emitted and the bath is empty. In terms of the system alone, when a trace over the bath is performed, this is 
a mixed state. Note that the mixed state is produced by ignoring the state of the bath. As far as future evolution 
of the system is concerned there is no classical uncertainty in the state of the system at this time, each of the ket 
vectors in the density matrix undergo a separate conditioned evolution, determined by the fact that the system must 
undergo one and only one emission during the interval \t\, t\ — T m ) (and also by any future measurement results). 

It is a surprising fact that despite these complications it is still possible to derive quantum trajectories for non- 
Markovian processes. 

In order to deal with this sort of mixed conditioned state we generalize the notation of Eq.(||) and write 

il RI (t,s)o\^s))=J2^Ri(t^W(s)), (12) 
k 

where 

n k RI (s) = (R(t,t )\U int (t,s)\k), (13) 
|V>( S )) = {k\U int (s,to)\I(t,t ))\i>(t )). (14) 

and the tilde denotes that the ket in unnormalized. Note that for any time s there is always one state of the bath, 
\I(t, s)), that is independent of the measurement results. Let us then take this state as one of the basis states of the 
bath, |fc = 0) = \I(t,s)). 

Having an explicit form for fl means that we are able to calculate the probabilities for the results of measurements 
made on the output in terms of the conditioned system evolution alone. We will now relate this to simulating the 
evolution of the measurement process forward in time. To do this in any practical way a necessary requirement is 
that there exists a time interval T m such that a measurement made at time t will be independent of the state of the 
system before time t — T m . In other words, there must exist a finite memory-time for the open system. We will discuss 



this in more detail in Sec. ID 



Consider a single run in which we have one particular input state and a record of the results of the measurements 
during the interval [t, to) and we would like to know the probability of a particular outcome during the next infinitesimal 
interval. The conditional probability density of getting a particular result during the next infinitesimal interval [t+dt, t) 
given the previous results over the interval [t, to) is 



E 



Piprfx^^ip^nmni P[R(t + dt,t Q )l{t + dt,t )} 

P[R(t + dt,t)\R(t,to)I{t,t )\= > 7X7771 m ' ( 15 ) 

j(£5t) P[R(t,h)I(t,to)\ 

P[R(t + dt,t )\I(t + dt,t )] 



,„_,„ P[R(t,to)\I(t,t )] 



P[I(t + dt,to)\I(t,t )]. (16) 
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If in a simulation we choose the input states consistent with P[I(t + dt, t)\I(t, to)] for each run then P[R(t + dt, to)\I(t + 
dt, to)]/P[R(t, to)\I(t, to)] will generate results R(t+dt, t) with the required probability distribution. These conditional 
probability densities can be written in terms of the operation as 

where we have defined the system state vector \i/j°(t + dt)), at time t + dt, as 

\iP°(t + dt)) = |^°(t)> + Ndn RI {t, t - T m ) o \^(t - T m )), (18) 



where \ip°(t)) = fl RI (t,t - T m ) o \^(t - T m )) and TV = 1/J (ip° (t)\ip° (t)) , again the tilde denotes that the state is 
unnormalized. This state determines the probabilities for the next measurement results. It is a transient state, in 
that it can be overwritten by the operation corresponding to subsequent measurement results (it is only the state 
of the system at time t if the measurement results from t to t + T m correspond to a state equivalent to the input 
state during this time). In Eq.([l8|) we have explicitly used the idea that the probabilities of getting a particular 
measurement result at time t do not depend on the state of the system or bath before time t — T m . If this is true then 
the system state at time t — T m is then p^y S (t — T m ) = J2 k i\(f — T m )\ip k (t — T m )) {ip k {t — T m )| and is independent of 

the measurement results at time t + dt. Here, Pk(t — T m ) = A4(ip k (t — T m )]ip k (t — T m )), where M. is a constant which 
ensures that Pk{t — T m ) = 1. This conditioned state converges to the reduced density matrix for the system when 
averaged over many runs of the above trajectories, p sys (t — T m ) = YIir /°sys(* ~ T m ). In the following sections we will 
investigate the nature of f2 for a particle conserving coupling and some quite general input states and measurement 
processes. 

B. System- bath interaction 

We assume that the interaction is linear and conserves particle number and also, for simplicity, that the bath is 
coupled to a single mode of the system. We consider situations that have bath and interaction Hamiltonians of the 
form 

tfbath = J dko; k 6 t (k)6(k) ) (19) 

ffi„t ^i^ljdk {. 9k 6t (k)a - ffk 6(k)at } , (20) 

where we have put Ti = 1 and a and 6(k) are the annihilation operators for the modes of the system and bath. The 
bath modes satisfy Bose commutation relations [6(k), 6^(k')] = 5 (k — k'). Note that our treatment does not depend 
on the system being a Bose field. The mode operator could equally well be replaced by, for example, the lowering 
operator of a two-level atom. The frequency of the kth mode is w k and is an arbitrary function of k. The effective 
coupling constant, g k , is normalized so that J dk|g k | 2 = 1, and the strength of the coupling is given by ^fy. The 
spatial dependence of the bath-system coupling determines g k , 

V^7k = J dr K (rK(r), (21) 

where /c(r) is the effective spatially dependent coupling constant and it k (r) are the spatial modes of the bath. An 
essential assumption in the following derivation is that the system and bath only interact over a finite region R 
localized about the origin so that we can put k(t) w for |r| > R. This is necessary so that we can later make the 
finite memory-time approximation. 

We can define input and output bath fields Jl5|,|3^] outside the interaction region (|r| > R) by 

6 in (i, r) = J dkb(k, t )e^ k(t - to) w k (r), (22) 

where to < t can be taken as the initial time and 

& out (i,r)= /" < flc&(k,ii)e-*' k < t - tl >Uk(r), (23) 
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where t\ > t can be taken as a time in the distant future. 

The unitary evolution operator Ui n t(t, to) is then given by the total Hamiltonian Eq.(Q) in the interaction picture, 



U in t(t,t ) = Texp I i^j / ds \^(s)ai(s) - £(s)a}(s) 



where 



£(t) = J dkg k b(k,t )e 



-iu> k (t—ta) 



(24) 



(25) 



is the driving field in the Langevin equations of motion for the system mode. a>i(t) is the system mode operator in 
the interaction picture and T denotes that the operators in the exponential are time-ordered. We define a memory 
function f m (t — t') as the commutation relation of the driving field with itself at an earlier time 



(26) 
(27) 



C. Measurement channels 



The quantum trajectory formalism presented here is based on an idealized situation where all the output from the 
bath is eventually measured by a perfectly absorbing measuring device. The physical situations that can be modeled 
by this formalism are not limited by this, as the information acquired from any measuring device that is not really 
there in the physical situation can be averaged over a number of runs of the simulated measurement process. So this 
formalism, although based on measurement theory, contains as a special case even the situation when no measurements 
are being made at all and the system is simply decaying into the bath. 

To describe separate measuring devices it is useful to introduce the concept of a complete set of measurement 
channels. A measurement channel has an associated field that contains only a range of modes of the original bath 
field. A measurement made on a particular channel is specific to the particular range of bath modes defined by the 
channel. The channels form a complete set, in that, together they contain the entire bath field. 

In a three-dimensional situation the bath output in all directions is being measured by detectors distributed around 
the system. To model the situation where each detector is sensitive only to the modes propagating towards it we 
can define channel fields containing only the appropriate bath modes for each detector. This concept of a channel 
also stretches to the case of spectral detection, i.e., measurements that directly access spectral information about the 
output field in a particular direction. A set of measurement channels can be defined such that each channel field 
contains only a range of frequencies of the output field. In this way, the results of a measurement made on one of 
these channel fields is specific to the particular range of frequencies defined by the channel. 

As a by-product, this formalism allows for the consideration of measurement devices that have a non-trivial spectral 
response. We can do this by choosing one channel to have the required frequency response, keeping its measurement 
record, and averaging the measurement results of the other channels (that make up the complete set) over a number 
of runs. 

In summary, we are considering a model where the output from the system passes through a generalized spectrom- 
eter, is scattered into different measurement channels and is then measured. The measurements are assumed to take 
place instantaneously at a single spatial point for each channel. The theory as it stands does not model the situation 
where the particles can scatter off the measuring device and reinteract with the system. 

The generalized linear-passive spectrometer is modeled by splitting the original output bath up into channels. A 
schematic of the situation is shown in Fig.^[ The number of input channels to the spectrometer must equal the number 
of output channels. 

We define the nth channel output field by, 

c° ut (^,r) = J dkS , n (k)6(k,ti)e- <Wk t*-*^tt k (r) + other baths. (28) 

where S n (k) is the spectral response of the channel. The contributions from other free baths (with the same dispersion 
relation as b(k)) do not interact with the system and are initially in a vacuum state. They do not concern us in latter 
derivations but they are necessary to preserve commutation relations (in other words they have a physical origin). 
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The definition of the channel fields in Eq. ( ]2q ) is in fact a linear transform of the mode operators of the independent 
baths, b m (k), to the mode operators of the channels, c n (k), for each k (we have written k as k for ease of notation). 
Writing bi(k) = b(k) for the original bath modes and defining c„(fc) by c™(t,r) = J dkc n {k)e~ lUk ^~ to ^ 'Wfe(r) , we can 
write the transformation as 

cn(k) = SraWhik) + J2 S nrn (k)b m (k). (29) 

m#l 

The requirement that the channel fields are independent and preserve commutation relations, [c„(fc), cj l ,(fc')] = 
5 n n'5(k — k'), is expressed as a constraint on S nm (k), 

^ Sjim{k)S mn r (fc) = 5 nn > (30) 

I a 

for all k. A special case of this when n — n' = 1 is 

J2\S m (k)\ 2 = l. (31) 

m 

where we have written S n (k) — S n i(k). This expresses the fact that the channel fields are a complete decomposition 
of the bath field. 

The nth output channel field is given by a unitary transform of the nth input channel field (the portion of the 
initial input bath that contributes to the measured output), 

c° ut (i,r) = Ul t (t,t )C(t,r)U int (t,t ), (32) 

= Vl/ dsh n (t-s,r)a{ S ) + C(t,r). (33) 



Note that the channel output fields contain a contribution from the system and the input channel fields. 
We have defined c'"(£,r) as, 

c™(t, r) = J dkS n (k)b(k, t^e-^^u^v) + other baths (34) 

and have also defined the impulse response functions, h n (t — s, r), for each channel, as the commutation relation of 
cjf (t, r) with the driving field £'(s) at an earlier time, 

h n (t-t',T)=[C(t,r),^(t% (35) 

dkS^k^ke-^-'VM (36) 



where t > t' . As the nth channel field is always measured at the same point in space, r„, we introduce the notation 
h n (t — s) = h n (t — s, r„). Definition Eq.(p^) yields the intuitively obvious relation with the output field (see Fig||), 

c° ut (t, r) = J dsh n (t - s)b out (s, r) + other baths. (37) 



The functions h n (r) and / m (''"), defined by the commutation relations Eq.(|35|) and Eq.(26), completely determine 
the characteristics of the system-bath interaction and measurement process. In the next subsection we will consider 
some general properties of these functions. 

D. Finite memory-time 

We now introduce an essential concept, that of the existence of a finite memory-time. That is, we assume that 
both the memory function and the impulse responses are such that at some finite time in the past T m the following 
conditions hold 
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t-T m 
t 

t-T™ 



dsf m (t ~ s)a(s) 
dsh n (t — s)a(s) 



> 



t—T m 



la 



t-T m 



dsf m (t - s)a(s) 



dsh n (t — s)a(s) 



(38) 
(39) 



for all n and all t during the entire measurement interval. To guarantee the above condition we assume that the 
memory and impulse functions satisfy 



ds\f m (t-s)\ 2 » / ds\f m (t-s)\ 2 , 

t-T m Jt 
t rt—T m 



ds\h n {t-s)\ z > 



ds \h n (t - s)| z 



(40) 
(41) 



which are independent of a(s). These are useful working conditions and are satisfied by all the functions considered in 
this paper. However, it is also possible to satisfy Eq.(|38|) and Eq. ( |39| ) if the memory function an impulse responses 
have a rapidly oscillating behavior. In this case the contribution to the integral after some time t — T m may average 
to zero. 

To understand the consequences of a finite memory-time we consider the physical effect of the driving field on the 
system. The Heisenberg equation of motion for a(t) is 



da(t) 
dt 



7 f* 

i[H sys , a(t)] - - fm(t - s)a(i 

1 Jt-T m 



(42) 



where we have neglected the smaller part of the integral. Correlations of the driving field with itself at earlier times 
cause the evolution of the system to be dependent on its past behavior. The existence of a memory-time implies that 
there is a finite time after which one can neglect the effect of the past state of the system in determining the present 
evolution. 

The physical effect of the system on the output channel field is given by 



t (i,r)«V7 



t-T™ 



dsh n (t - s, r)a(s) + cj"(t, r), 



(43) 



where we have again explicitly inserted the lower bound of the integral. From this equation one can see that the state 
of the system in the past determines the output channel field in the present and that the finite memory time imposes a 
cutoff time for this dependence. This relation highlights an important point regarding the interpretation of the results 
of measurements of the output channels in terms of the system dynamics. A measurement result at time t does not 
mean that we can ascribe a state corresponding to this result to the system at the same time t. Unmeasured output 
at later times is still quantum mechanically entangled with the system at time t and a subsequent measurement could 
require us to revise our knowledge of the system state completely. In the Markov case, because of the delta function 
impulse response, a measurement result meant that we could specify the state of the system at the same time. In 
the general case considered here, only states of the system at a time t — T m can be thought of as fully specified, as 
the output after time t is uncorrelated with the system before this time and future measurements will not alter our 
knowledge of the system state. We can then interpret the state at time t — T m as the system state conditioned on the 
measurements up to and including time t. 



E. Conditioned system evolution 

We now return to a treatment in terms of state vectors and set about determining the explicit form of £Iri for a 
general measurement scheme. 

We consider quantum trajectories to be an attempt to simulate the measurement results of a distant measuring 
device by the evolution of a local system emitting particles. This is related to the way in which the spectrometer is 
treated in this work. Instead of explicitly modeling the generalized spectrometer as a physical system (coupled to our 
localized system) we move towards a signal processing treatment, where the spectrometer is treated as a black box 
and its inputs are related to its outputs via a spectral response. 

Take, for example, a number of measuring devices making particle counting measurements. At any one time we 
assume that there are only two possible outcomes (for perfect detectors): the null result, where no particle is present 
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at any of the detectors; and the positive result, where a particle is detected at a single detector. A positive result 
corresponds, in general, to a superposition of particle emissions from the system at times in the past. Therefore, a 
probability can be calculated for the event from only the evolution of the system undergoing emissions. We assume that 
the weighting given to each emission time characterizes completely the passage of the particle to the detectors. The 
probability of the null event can be inferred from the requirement that the probabilities of all possible measurement 
results must add up to one. The probability of the null event cannot be determined directly as we are not explicitly 
modeling the physical evolution of the particle from the system to the detector. 

It is a rather interesting fact that in the spectral detection case the null result, the no detection result, cannot 
be assigned an operation. In the Markovian quantum trajectory case a null result corresponds precisely to a non- 
emission event (a virtual emission and absorption) and so an operation is also assigned to the null event. The virtual 
emission-absorption event is localized in the interaction region and so is not directly related to results at the distant 
detector. However, in the Markovian case the detector can in principle be brought right up against the interaction 
region (which takes place at a point in the Markovian case). 

Recall that we have defined the operation in the system space corresponding to the measurement results over the 
interval [t, to) as 

n«/(*,*o) = (R{t,t )\U ta t(t,to)\I(t,t )), (44) 

where U- m t(t, to) is given by Eq. (|24|) . We are trying to simulate the measurement process in terms of the evolution of 
the system alone and so have not included the spectrometer in the unitary evolution. The operation contains all the 
information about how to evolve the system forward conditioned on measurement results and also how to calculate 
the probability of a particular outcome. 

We assume that our idealized measuring devices yield a result for each infinitesimal interval of the continuous 
measurement process. The state corresponding to the results must then factorize into a product of states over each 
infinitesimal interval. We can therefore write the output state as a product of state generators over each infinitesimal 
interval acting on the vacuum. We consider measurements where the generator for the output is written as a function of 
the creation operators for the free input channel fields at a particular position c^ n (t, r„), (the measurement apparatus 
can be thought of as a sink for particles at position r„). In mathematical language we have, 

\R(t,to)) = n G M n fe-rn)]|{0}>, (45) 

where Gf n [c^ n (tj, r„)] is the generator of the state during the interval [tj + dt, tj) where tj ranges from to to t — dt. 
The input channel operator appears in this expression because it is the interaction picture output channel field. 
For example, these generators may correspond to measurements of the quadrature amplitudes (homodyne O and 



heterodyne detection \22\) or absorption measurements (particle counting) |2||. 

The input can similarly be thought of as a source of particles at positions q and is written in terms of the input 
field, i.e., 

\I(t,t o ))=l[G I J [b^(t 1 ,c l )}\{0}). (46) 

3 

For example, these generators can correspond to vacuum or coherent input states, squeezed input states or thermal 
input states 

Recall that our problem is to determine the operation Q,Rj(t,to) in terms of the system variables only. Because of 
the factorized form of the input and output state we can consider operations of the form 

O n (t 3 ) = ({0}|Gf„[ C ^(t J ,r„)]C/ int (t,io)G J / [6t i "(^,q)]|{0}), (47) 

for each tj and each channel n. 

To proceed further we note that we can write the interaction unitary evolution operator t/i n t (£, io) m the form 

U bA (t,t )=Ta{V+VoV-}, (48) 

where T a denotes time-ordering of the system operators a/ and a\ and where 

y + =exp jv7^W(sK(s)}, (49) 
Vq = expj-~^ dS2 J t dsifmfa - si)<4(s 2 )a/(si) j , (50) 
V- =cxp(-V7 f dst;( S )a\( S )} . (51) 
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This can be proved by differentiating Ui n t(t,to) and re-ordering the bath operators. The operation O n (tj) then 
becomes 

0»(*i) = ({0}\T a {GfJCit^TnW+VoV-G^itj,^)}} |{0}), (52) 

where we have expanded the brackets of the time-ordering to encompass the state generators as they are only a 
function of the bath operators. 

We now make use of the relation (compare with Eq.(|2"s|)) 

F_&t-(t, q )y_-i =b^(t,q)+^(t,t ;s), (53) 
V^C(t, v n )V+ = C(t, r„) + ^An(t, t ; s) (54) 

where we have defined for notational purposes 

A n (t, t ; s) 
A(t, t ;s) 

We then make the above replacements and note that V_ can be annihilated on the right-hand vacuum and V+ on 
the left-hand vacuum. The the system operators aj and a\ commute with all the input fields c™(tj, r n ) and b ln {tj 7 q), 
as they are still inside the time-ordering brackets and so act like complex numbers. 

If we repeat this procedure for each time tj and if we restrict the nature of the state generators so that when the R 
and / generators for each time are swapped from left to right no bath modes are left that are not normally-ordered 
(specific examples will be given latter) then we are left with an operation of the form 

fi«z (*,*>) = rj YlGf^Anitj^sWollG^A^tj^s)} \ , (57) 

where we have dropped the reference to the mode operators in the time ordering operator T . 

This is the general result for arbitrary (with the above restriction) measurement schemes and arbitrary input. No 
assumptions have been made regarding the strength of the coupling constant or the nature of the correlations between 
the system and the bath after the initial time. This operation can be transformed to the Schrodinger picture by 
multiplying on the left by the unitary evolution operator for the free system. 

Because of the complex nature of the above operation the measurement process can only be simulated on 
a computer in the weak-interaction short-memory-time case where I(t)/T < 1, for all times of interest, where 
I(t) = (6 tout (i)& out (t)) ~ 7(a t (^)a(i)) is the intensity of the output. T ~ 1/T m is some measure of the rate of 
decay of the of the memory function or impulse response. It determines, in essence, how non-Markovian the process 
is. In this situation we can make a perturbative expansion of the operation f2flj in powers of "f(a' (t)a(t)) /T. To lowest 
orders the dimensionality of the integrals during a memory-time will be small and they can be evaluated numerically 
with a computer. In absorption measurements each click of the counter corresponds to one integral so this can be 
simulated in the case where there are on average only a few detections (real emissions) during a memory-time (there 
will also be integrals of the same dimension corresponding to virtual emission- absorption events). We shall explicitly 
demonstrate how this can be done in the case of optical spectral detection in the next section. 

II. OPTICAL SPECTRAL DETECTION 

In this section we consider a special case of the above theory; that of spectral detection of an optical field emitted 
from a source. To check the accuracy of our numerical simulations we shall apply it to a known problem; the spectral 
measurement of light emitted from a strongly driven two-level atom. 

Spectral detection of an optical field is defined here as meaning that each photodetection can be identified with the 
emission of a photon from a certain part of the source spectrum |]l9| (this is in contrast to measurements from which 
the spectrum can be computed, such as heterodyne detection). Spectral detection necessarily involves the interference 
of light emerging from the system over some time period At. The length of this time period varies inversely with the 
accuracy of the spectral measurement. 



dsh n (t — s)ai(s), 



to 



d S [b^(t,q),^( S )}ai( S ). 



(55) 
(56) 
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The type of spectral detection that has been considered previous to this work is that of detecting fields with 
spectral peaks of widely differing frequencies. These are then just treated as independent Markovian baths for the 
source to decay into. The evolution of the wave function for this type of system can then be treated with the standard 
quantum trajectory approach. Here we would like to go beyond this and treat wave function evolution for the source 
conditioned on spectral detection of emitted photons within an emission linewidth. To achieve this it is necessary 
to give up knowledge of the exact time of emission. We can only know the time of emission to within an interval 
At > 1/Acj if we are to resolve the emission spectrum to within Auj of a certain frequency. 



A. The spectrum of a strongly driven two level atom 

A strongly driven two-level atom can be described (in a frame rotating at the frequency of the driving field) by the 
free system Hamiltonian, 

i? a tom = ^(<T t +*), (58) 

where a is the lowering operator and f2 is the classical driving field strength and we have assumed that the driving 
field is on resonance with the atomic transition. 

The optical field interacts with the atom in a very small region about r = and we can write the unitary evolution 
Eq.(||) as 

EWMo) = Texp [ { VlJ t ds [*■}(*)£(») - <7/(s)£0)] | , (59) 

where in this case 

POO 

£(*)=/ do&(cj,i )e- ,iw(t ~* o) . (60) 



We have replaced the mode label k by the angular frequency cu as is customary in optics in the situation when 
the direction of the emitted light is unimportant. We have also ignored the polarization of the light field. Note 
that £(t) is the same as the free input field at the point of interaction r = and it is delta-correlated with itself, 
C (*')] = — O- Because of this delta function correlation of the driving field a Markovian master equation 
can be derived for the system by tracing over the bath (sec |jlf for example). Here we will consider a non-Markovian 
unraveling M of this master equation into quantum trajectories . 

It has been shown by Mollow p3] that the spectrum of the spontaneous emission from a strongly driven two-level 
atom has three peaks at the frequencies 0, ±0. It has also been shown by Cohen- Tannoudji and Reynaud [^0j that 
the side-bands exhibit photon antibunching independently, but are correlated via photon bunching. Experiments by 
Aspect et al. [E7j have confirmed these theoretical results by frequency filtering the light incident on three photodetec- 
tors, so that each photodetector is only sensitive to a range of frequencies about one of the peaks of the output light. 
Here we would like to reproduce these well known results by simulating this "spectral detection" of the spontaneously 
emitted light with a non-Markovian quantum trajectory for the source atom. 

The original theoretical predictions were made with a fully quantized atom-field model and the model was analyzed 
in terms of the eigenstates of the Hamiltonian, the dressed states of the atom. Instead, we will consider the simplified 
model described by the Hamiltonian above, Eq.(|5^). In the limit of a large classical driving the field can be assumed 
to be in a coherent state and our model is valid. All the important results of the dressed state model will be present 
in our treatment Jl9| . We can qualitatively understand the behavior of this system by considering the evolution of 
the interaction-picture lowering operator, 

<rj(t) = + l+)(-|e in * + |-><+|e- ia *) (61) 

where |+), |— ) are the eigenstates of a x = <r+ + er_, cr x |±) = ±|±). If we are continuously monitoring the emission 
from the system and are able to resolve the frequencies of the emitted photons, then if a photon is detected in the 
+fi(— fi) frequency peak the system will be projected into the opposite eigenstatc |— )(|+)). If, however, a photon is 
detected in the central or zero frequency peak the value of (<J x (t)) will be unchanged. The interesting behavior of this 
system stems from this simple idea. 

In the following sections we will only consider the case of direct detection, that is, absorption measurements by 
photodetection. The state corresponding to an absorption measurement in one of the spectral channels during an 
infinitesimal interval is 
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\R(t + dt,t)}) =l[c^(t) dN ^\{0}), (62) 

a 

where the restrictions (J2 n dN n (t)) 2 = ^2 n dN n (t) and dN n (t) 2 = dN n (t) describe a point process, such that in an 
infinitesimal interval, cither only one photon is detected in one of the channels or no photons are detected in any of 
the channels. The dN n (t) can take the value 1 or corresponding to a click of the photodetector for channel n and no 
click of the photodetector. We also assume that the input state is the vacuum (the coherent driving being included 
in the system Hamiltonian). 

Consider the situation where N photons have been detected at each of the times t\ < t 2 < ■ ■ • < U < ■ ■ ■ < tjv m the 
channels m, n.2, ■ ■ ■ , rii, . . . , njv during the past interval [t, t — T m ). The wavefunction Eq.(^)for the next infinitesimal 
interval is 

t 7 



\ft(t + db)) = \^(t))~dtla\(t)ajm (t)) 



^dN n (t) 



2 r t 



^VlY] / dsh n (t ~ s)T {*i(8)A% t - T m )} \tf{t - T m )> - |^°(t)> 



(63) 



where 



and 



A*(t, t - T m ) = Y[ s/lAn, {tj,t - T m ; Sj)V (t, t - T m ), (64) 



ft j 

A nj (t 3 ,t-T m ;s) = / dsh nj {tj - s)a I (s). (65) 

Jt-T m 

In this case Vo(t,t — T m ) = cxp{— ^ f t _ T dsa\(s)ai(s)}, this leads to a much simpler evolution than the general case. 
The integrals A nj (tj, t — T m ; s) have the lower bound t — T m . However, there is no factorization of the operation at 
t — T m . Instead, we must write a different state vector for each A 1 . The product is over each combination of from to 
N integrals. To understand what we have done here it is necessary to look at the integrals more carefully. Consider, 
as an example, two time-ordered integrals, 

T{A ni (ti,ti-T ni , Sl )A n2 (t 2 ,t 2 -T m ,s 2 )V (t,t2-T m )}= f % [' T {■••}, (66) 

where t > t\,t% > t — T m . For each detection at time U there is an integration over emission times during the interval 
[ti,ti — T m ). If we split both integrals into the sum of two parts - emissions before t — T m and emissions after - we 
can write this as the sum of four terms 

T{-..}= / T{..-}T{...}+ / T{.-.} T{-..} 

t 2 -T m Jt-T m Jt-T m Jt-T n , Jt 2 -T n , 

T{...}/ T {..-} + T{...} / / T{-..}. (67) 

t-T m Jt!-T m Jti-T m Jt 2 -T m 

This can be easily generalized to N integrals. This is the essence of what has been done in Eq.(|63|). The parts of 
the integrals before t — T m get absorbed into wave functions. Each wavefunction is labeled so that its subsequent 
evolution is consistent with its previous evolution. Physically this means that if a photon is detected in a channel at a 
particular time this corresponds to a sum over individual emission events and a photon cannot be emitted more than 
once. This is the detailed structure glossed over in the o-product notation of the previous section. 

The probability density for a detection to occur in channel n during the time interval [t + dt, t) conditioned on the 
previous detections is then given by 

P n (t + db\tif,t N - U ~.,ti) = + (* + *)). (68) 

where 



2 N i 

\i%(t + dt))=ffSyJ2 f dsh n (t-s)T{aj(s)A i (t i ,t i -T m )}\i; i (t-T m )). 



(69) 
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This result is of limited practical use in the above form as evaluating multi-dimensional integrals on a computer is a 
slow business. Instead we take a less ambitious approach and assume that the bath-system coupling is weak so that we 
can make what amounts to a perturbative expansion of the trajectory in terms of the coupling strength y/j. In other 
words, we assume that on average only a small number of detections occur during the time interval [t, t — T m ), the 
number of detections determining the dimension of the integration. To first order, i.e., when the maximum number 
of detections per memory time is one, the trajectory is very similar to the Markovian case. We therefore consider the 
next highest order so that there is a finite possibility of two detections occurring per memory time. 

In order to demonstrate how a trajectory is simulated on a computer in practice, we write down as an example 
the Schrodinger picture form of Eq. (|6^) in the case where one detection occurred at time t' in channel n' during the 
previous memory time, 

' ' '"^ = 7 / dsi [ ds 2 h n (t-s 1 )h nt (t'-S2)U eS (t,s 1 )aU eS (s 1 ,s 2 )aU eS (s2,t-T m )\i;%t-T m )) 



M Jf Jt-T a 

+"fU cS (t,t') I dsi I ds 2 [h n (t~ s^h^it' - s 2 )} synl U c g(t\s 1 )aU cS (s 1 ,S2)aU cS (s2,t-T m )\%jj (t-T m )) 

Jt-Tm Jt-T m 

+V7 / dsh n {t - s)U cS (t, s)aU eff {s, t - T m )\^\t - T m )>. (70) 

Jt-T m 

where U eS (t,t ) = exp{[-iff atom - a](t - t )} and [h n (t - s\)h n >(t' - s 2 )] sym = h n (t - si)h n /{t' - s 2 ) + h n (t - 
S2)h n '{t' — si). The right-hand side is the sum of three terms. The first term represents the situation where the 
emission corresponding to the detection at t' occurs after time t — T m but before the emission that is detected at 
time t. The second term represents the situation where both the t emission and the t' emission occur during the 
same time interval [t',t — T m ) and the time-ordering of the emission is accounted for by the symmetric ordering of the 
impulse response functions. In the last term the time t' emission has occurred before time t — T m (and is included in 
\^} 1 {t — T m ))) and only the time t emission is left. The states \ip l (t — T m )) are propagated forward one time step dt 
by doing the integrations contained in the operation from t — T m to t — T m + dt. 

An important technical aspect of the simulation is; how to choose a finite length impulse response for numerical 
calculations given a particular spectral response? This is a problem of signal analysis and numerical routines exist 
for automatically generating finite length impulse responses from a given spectral response such that the spectral 
response corresponding to the finite length impulse response is as close as possible to the original in a least-squares 
min-max sense ||J. Note that if the impulse response is too long there is a greater chance of having to deal with 
higher-order integrals, while if it is too short the spectral response will not approximate the original spectral response 
and the probabilities will not be accurate. In line with this we have avoided triple integrals by shortening the 
response functions when calculating the probabilities for detection after two detections have already occurred during 
the previous memory-time. As the probability of three detections during a memory time is already small (in the limit 
of weak damping) this should not have a large effect on the statistics of the counts. 

The simulation procedure at each time step dt is as follows. First, the conditional probabilities for detections at 
time t + dt in each of the channels during the time step is generated by calculating the wave function Eq. (|7^) for each 
n. This depends on any previous detections during the time interval. The normalization J\f = (ip° (t)\i(} (t)} must also 
be calculated independently at each time step by propagating from t — T m to t conditioned on previous detections but 
assuming that there will be no more detections during the time interval [t + T m , t). The probabilities for detections 
in each channel are added together and compared to a random number. If the probability of a count is less than the 
random number, the state at time t — T m , | ip^ (t — T m )) is propagated forward conditioned on the fact that there was 
no detection at time t + dt. If the probability of a count was greater than the random number a channel is randomly 
picked from the probability distribution over the channels. The state at time t — T m is then propagated forward dt 
conditioned on this outcome. The process can then begin again. 



1. The Frequency Filter 

Here we consider the situation where we frequency filter the output from the atomic source. The light emitted from 
the source passes through a Faraday isolator (so the filter cannot affect the source) and then through the filter. It is 
therefore possible, by tracing over both the bath and the filter cavity, to write down a master equation for the source 
alone. 

We consider as a frequency filter a two-sided Fabry-Perot cavity that can reflect and also transmit light. The 
coupling strength n is assumed to be the same at each mirror. The transmitted light will be in a Lorentzian shaped 
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band of width k about the resonance frequency of the cavity, v. This light is detected by a photodetector as is the 
light reflected from the other mirror. 

The reason that the non-Markovian quantum trajectory for the evolution of the atom can generate the correct 
probabilities for the spectral detection without explicitly modeling the filter cavity (with a Hamiltonian) is because it 
is possible to eliminate the filter cavity mode completely and write the output from the cavity in terms of the input 



We define two output channels, c™ and Cj . , as the two output fields from the reflecting and transmitting mirrors, 
respectively, after the output from the atom 6 out has interacted with the filter cavity, 

4 ut W = 4(Mo)[& out W + c^(t)]U F (t,t ), (71) 

= r t (t) + c! r n (i) + ^a(i), (72) 

Cr ut W = U F (t,t )[b out (t) + c?(t)}U F (t,t ), (73) 

= b out (t) + + vM*), (74) 

where Uf is the unitary evolution for the interaction of the baths with the filter defined by 

U F (t,t )=Texp^i^K^ ds (a 7 (s)[4i n (s) + fe tout (s)] - H.c?j J , (75) 

where aj(s) is the cavity mode in the interaction picture, ctp n (s) is the free transmitted field and b° ut (s) is the output 
field that is coupled to the atom, both are delta correlated. 

We can now easily solve the Heisenberg equations of motion for the filter cavity mode by Fourier transforms. If we 
assume the state of the filter in the distant future is the vacuum state we can put a(t±) = 0, and we find 

c out (t) = i r d ^-^t,)+ K c T ^h) e _ M ^ tiK 

V27TJ-00 K-l\U-V) 

4 «t (f) = i r ^-vWjujt^+^M)^^ (77) 

These are the channel field operators corresponding to photodetections of the reflected and transmitted fields of the 
filter. The presence of the bath ct preserves commutation relations. The impulse response for each of these channels 
are therefore, 

4 ltcr (i -*) = 4= r du> %{L °7 V) - e~^-») , (78) 
V2W-00 k-%{u-v) 

= -5(t-s) + h%i tcr {t-s), (79) 
4 lter (t - s ) = -== / dio p- K e-^- s \ (80) 

= u(t - s)Ke- lu(t ' s) e- K{t - s) , (81) 

where u(t) is the unit step function, zero for r < 0. Note that in line with the general discussion in the introduction 
if we detect a photon in the channel c!^ ut then we know its frequency to within the cavity linewidth k but have lost 
information about the time of emission to within \j k. 

We can simulate the spectral measurement process via Eq.([33|) by simply substituting in the filter impulse responses. 
The perturbation expansion in this case is in orders of 7/K as k defines the decay time of the impulse response functions. 

Fig.|| is a schematic of the situation in frequency space, showing the three peaked output and a superimposed 
cavity filter centered on the central peak. In this situation the central peak is transmitted while the two side peaks 
are reflected by the filter. 

The results of a simulation with the parameters of Fig.|| are shown in Fig.|| where we have plotted the photodetection 
waiting times for the reflected and transmitted emission for one run of a total of 12 x 10 4 detections. The time average 
of a single realization is equivalent to an ensemble average. The transmitted light is mostly from the central peak and 
the reflected from the side peaks. Notice that both these waiting-time distributions show a marked increase in the 
frequency of longer waiting times compared to the distribution for the combined emission. As 7 ~ k the resolution 
of the spectral detection was of the order of the decay linewidth. This is reflected in the fact that a non-vanishing 
fraction of the photons were detected within a memory time T m = I/7 of each other. 
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A plot of the evolution of the probabilities of the two channels for the first ten detections of a single trajectory is 
shown in Fig.||. The probability of a transmission is proportional to the expectation value of the number of photons in 
the filter cavity. When there are no detections it oscillates at the frequency O, in time with the oscillation of (a y (t)). 
The probability of detecting a reflected photon is an interference between the possibility of a photon being reflected 
directly off the mirror and the possibility of a photon coming back out of the cavity through the reflecting mirror. 
The oscillations are suppressed by this interference. The possibility of detecting a reflected photon is very sensitive 
to the phase between the states representing these two possibilities, an example is the large peaks in the reflected 
probability when a photon is detected in the reflected channel and the oscillation has the right phase. 

Fig.^ is a plot of the time evolution of the conditioned evolution of the expectation value of a x and a y during the 
same trajectory. Note that the conditioned system shows signs of having started to emit even before a detection occurs. 
This is clearest in the decaying (a y (t)) oscillation before the second and sixth detection. This is a clear indication of 
the non-Markovian behavior of the conditioned state; measurement probabilities depend on the previous state of the 
system and the conditioned state of the system is affected by a detection at a later time. The interpretation of this 
behavior goes as follows: at the time of detection the system + bath is in a superposition of all the possible states, a 
detection selects out the state that corresponds to one particle in the bath; associated with this state is a particular 
history, i.e., a weighted sum of emission times; this history then determines the conditioned state. In the filter cavity 
case a x is a constant of the free motion and is unaffected by the measurement process as the measurement does not 
distinguish between the sidebands. 

This situation can also be written in terms of the theory of cascaded systems |2^,^t|. In that case, the evolution of 
the filter cavity as well as the atomic system are treated together as a coupled system. This is a special case of the 
more general idea that was mentioned in the introduction. The atom and the filter are coupled to the reflected bath 
at the same physical location and they are also coupled together directly. The quantum trajectory is then evolved 
by simulating the conditioned evolution of the complete system of source and filter. This cascaded system trajectory 
can be described by the effective Hamiltonian (in the Schrodinger picture), 

H e g = -ff a tom + {v - iK,)a)a - i + V7^ crat ) i ( 82 ) 

with the collapse operators 

C T = v^a, (83) 

C R = v^a + yfytr. (84) 

Because we can formulate this situation in two different ways (in terms of normal quantum trajectories for an extended 
system, and in terms of a non-Markovian trajectory for the source alone), comparing the results of the two simulations 
provides a good check of the non-Markovian trajectory. In fact if one was interested in this particular system simulating 
the evolution of the extended system is much simpler and requires much less computer time than the non-Markovian 
equivalent. 

In Fig^ we have plotted the various waiting times taken from simulations of both the non-Markovian method and 
the Markovian method for the extended-system the results show good agreement, with the discrepancy within the 
statistical fluctuation. In general, the non-Markovian simulation takes of the order of 10 times as long to run on a 
computer. 

In the next section we will give an example where it is not obvious how to build an extended system that would 
accurately simulate the measurement situation. 

2. The Prism 

We consider a simple model of spectral detection performed by a prism, where the light emitted from a source 
propagates through a prism (or a spectral grating) which spreads the light into a spectrum and is then incident on 
an array of photodetectors. Each detector is then effectively sensitive to a sharply defined band of frequencies. We 
can then model this situation by assigning a top-hat frequency response function to each photodetector (labeled by 
the variable n) centered about a frequency ui n . The output channel field is given by 

-, /-w„+A/2 

c° ut (i) = -= \ duh{^M)er iu{t ~ tx) + other baths, (85) 

where A is the width of the band (which we have assumed are the same for all n). The prism then has the impulse 
response functions 



15 



(86) 



These response functions do not obey causality. This is because we have not included the propagation time from the 
system to the detector. A finite approximation (which is physically valid as an infinitely sharp cut in frequency is 
non-physical) to this will be nonzero over a time interval of [T m /2, — T m /2). We can then take the greatest time of this 
interval as the time of a detection. The prism example has the advantage that we can easily model the detection of 
the three peaks of the Mollow spectrum by assigning a channel field to each of the peaks. We consider a spectrometer 
that splits the light into three frequency bands with each band centered on a separate peak of the Mollow spectrum. 
A schematic of the situation is shown in Figj7[ 

The results of a computer simulation of 8 x 10 4 detections are shown in Fig.|], with system parameters as in Fig.|^ 
and with T m = I/7. The results show a definite anti-bunching of the side-peak photons (characterized by the peak of 
the distribution being shifted to longer waiting times). They also demonstrate that the inter-sideband waiting times 
(i.e., the time between an emission into one sideband and an emission in the other) have a distribution similar to that 
of the central peak. Note that this plot shows details of the waiting-time distributions closer to zero than the filter 
case. 

A plot of the detection probabilities in each channel for one run of ten detections is shown in Fig.|9|. We have 
labeled the channels; L (left), R (right) and C (center), referring to the side and center peaks of the Mollow spectrum 
that the channels correspond to. For our purposes the thing to note about these trajectories is that after a detection 
occurs in the right (left) channel the probability to get another detection in the same channel is nearly zero until a 
detection occurs in the left (right) channel regardless of detections that occur in the center channel. This is just the 
predicted anti-bunching in the side-peaks. This phenomenon is reflected in the conditioned state of the atom. A plot 
of the expectation values of a x and a y for this conditioned state are shown in Fig.[l0| Whereas in the cavity filter case 
(&x(t)) was zero throughout the trajectory, here, because the measurement process can distinguish the side-peaks, 
the atom gets projected into the eigenstates of a x with each measurement of a side-peak photon. The amplitude of 
oscillations in (a y (t)) decrease when \(cr y (t))\ is large as the atoms behavior is basically restricted to the surface of 
the Bloch sphere. Note that the conditioned expectation values (cr y (t)) and (cr x (t)) show even more obvious signs 
of future detections than the filter cavity case, for example, the very obvious decay of the oscillations of (cr y (t)) and 
the motion of (a x (t)} towards zero long before a detection. This is more pronounced in the prism case because the 
impulse response is a sine function as opposed to a more rapidly decaying exponential in the filter case. 



CONCLUSION 



We have derived a general form for non-Markovian quantum trajectories corresponding to a particle conserving 
coupling between the system and bath. We have also introduced the concept of a complete set of measurement 
channels to model measurement devices that are sensitive to only a range of modes of the bath. In the limit of 
weak coupling the non-Markovian quantum trajectories so derived are amenable to computer simulations. We have 
demonstrated the practicality of the method by making computer simulations of the special case of optical spectral 
detection. The results agree with the predictions of more traditional Markovian methods. The theory is general 
enough to be applied to a number of areas in physics where open systems arise and where the Markovian assumption 
cannot be made. In further work we hope to apply the above theory to simulating the statistics of the output coupled 
atoms from a Bose-Einstein condensate |l^,[l^|. Another topical problem to which these methods could be applied is 
that of radiation into band-gap materials ||12U13|1 . 
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FIGURES 




FIG. 1. Diagram of the weighting distribution over the times of emission for a particle that may be detected in the output at time 
t (dashed line) and that has been detected at time t\ (solid line). The times at which previous emissions can be neglected t — T m and 
tl — T m are also shown. 



INPUT 




FIG. 2. Schematic of the spectral detection setup. The output from the system of interest is an input to one of the spectrometer 
channels. The inputs to the other channels are simply the vacuum. The spectrometer mixes the channel inputs so that the output 
channels contain certain spectral components of the original system output (as well as contributions from the other channel inputs). 
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FIG. 3. Schematic of the filtering situation considered. The three peaked spectrum of the two level atom is shown and superimposed 
on top is the Lorentzian lineshape of the Fabry Perot filter. The filter cavity is on resonance with the atomic transition v = (the central 
peak of the Mollow spectrum) and has a line width larger than the atomic decay rate, k = 57. The side peaks of the Mollow spectrum 
are at ±£7, where Q = IO7. 
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FIG. 4. Plot of the results of a simulation of the filtering situation depicted in Fig.g, all three plots are waiting time distributions (in 
arbitrary units) determined by bining the times between each subsequent detection in 100 bins. The solid line is from the non-Markovian 
quantum trajectory, the dashed line is the Markovian method for the extended system, (a) is a plot of the combined photodetection 
waiting time distribution from both the reflected and transmitted emission, (b) shows the waiting time distribution for just the reflected 
emission (which consists mostly of the two side peaks), (c) is a plot of the waiting time distribution for the transmitted emission (mostly 
from the central peak). The ratio of the transmitted and reflected emission is also given. 
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FIG. 5. Plot of the probabilities for detecting a transmitted Pt and reflected Pr photon during the evolution of a single quantum 
trajectory for the first ten detections. The labels R for reflected and T for transmitted at the top of the graph denote the channel that 
the detection occurred in and the dotted vertical lines indicate the time of the detections. 
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FIG. 6. Plot of the conditioned expectation values for a y (solid line) and a x (dashed line) for the same trajectory as in FigJ^. 
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Prism 




FIG. 7. Schematic of the prism-like spectral detection situation. The three peaked spectrum of the two level atom is shown and 
superimposed on top are the top-hat functions corresponding to the detectors centered on each peak. The spread of the frequency band 
of each channel is A/2 = 57 and lu^ = — fij = +^ an d = 0, where, as before, f2 = IO7. 
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FIG. 8. Plot of the results of a simulation of the prism depicted in Fig.jj], all three plots are waiting time distributions determined by 
bining the times between each subsequent detection in 100 bins, (a) is a plot of the combined photodetection waiting time distribution 
from all three spectral bands, (compare with FigQ (a)), (b) shows the waiting time distributions within each spectral band (each band 
coinciding with one peak of the emission). The solid line corresponds to the central peak and the dashed and dotted lines to the side-peaks, 
(c) is a plot of the intcr-sideband waiting time distribution, i.e., the time taken between a detection in one of the side peaks and a detection 
in the other. The ratio of the emission into each peak is also shown 
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FIG. 9. Plot of the probabilities for detecting a photon from the central peak Pq and the two side peaks Pj, and Pjj during the 
evolution of a single quantum trajectory for the first ten detections. The labels C for central and L and R for the side-peaks at the top 
of the graph denote the channel that the detection occurred in and the dotted vertical lines indicate the time of the detections. 
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FIG. 10. Plot of the conditioned expectation values of a y (solid line) and a x (dashed line) for the same trajectory as depicted in Fig. 
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